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Abstract 

A highly accurate self-consistent particle code to simulate the beam-beam col- 
lision in e^e~ storage rings has been developed. It adopts a method of solving 
the Poisson equation with an open boundary. The method consists of two steps: 
assigning the potential on a finite boundary using the Green's function, and then 
solving the potential inside the boundary with a fast Poisson solver. Since the 
solution of the Poisson's equation is unique, our solution is exactly the same as 
the one obtained by simply using the Green's function. The method allows us to 
select much smaller region of mesh and therefore increase the resolution of the 
solver. The better resolution makes more accurate the calculation of the dynam- 
ics in the core of the beams. The luminosity simulated with this method agrees 
quantitatively with the measurement for the PEP-II B-factory ring in the linear 
and nonlinear beam current regimes, demonstrating its predictive capability in 
detail. 
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1 Introduction 



The beam-beam interaction is one of the most important hmiting factors determining the 
luminosity of storage colhders. It has been studied extensively by theoretical analysis 
experimental measurements 0, and computer simulations 0. Historically, due to the com- 
plexity of the interaction, many approximations, such as strong- weak ^ or soft-Gaussian ^ , 
have been introduced in order to simulate the interaction in a reasonable computing time. 
The self-consistent simulation of the beam-beam interaction by solving the Poisson equation 
with a boundary condition has been proposed first to investigate the round beams and 
then the flat beams [0. To enhance the accuracy and to reduce the computational over- 
head, an algorithm (and a code) of the so-called 6f method that can handle strong-strong 
interactions has been introduced Another self-consistent approach to the beam-beam 
interaction is to use the Green's function directly or indirectly p!0| . 

In the present paper we will develop a method that takes advantage from both self- 
consistent approaches: a smaller region of mesh from the method of using the Green's func- 
tion and a faster solver for the interior. In order to develop a highly accurate predictive code 
at the luminosity saturation region, it is necessary to have a fully self-consistent treatment of 
field-particle interaction at collision. Since we are interested in simulating the Asymmetric 
e~^e~ Storage Collider PEP-II [0, which needs to maximize the luminosity and thus the 
beam current, it is even more crucial that the beam-beam interaction in the large current 
regime be treated accurately. 

In a self-consistent simulation of the beam-beam interaction in storage rings, the beam 
distributions have to be evolved dynamically during collision with the opposing beam to- 
gether with the propagation in the rings. During collision, the beam distributions are used 
at each time sequence to compute the force that acts on the opposing beam. 

Since positrons and electrons are ultra-relativistic particles in high energy storage rings, 
the beam-beam force is transverse and acts only on the opposing beam. Hence, given a 
beam distribution, we can divide the distribution longitudinally into several slices and then 
solve for the two-dimensional force for each slice. Self-consistency is achieved by introducing 
many-body particles in the field that in turn constitute charge-current, the strategy of the 
particle- in-cell (PIC) procedure (for example, Ref. |T^). In this paper, for simplicity, we use 
only a single longitudinal slice for a bunch, ignoring any beam-beam effects encompassing 
over the length of the bunch. 

2 Method 

In modern colliders, beams are focused strongly at the interaction point to achieve high 
luminosity. As a result the transverse dimension of the beam is much smaller than the 
dimension of the beam pipe at the collision point. Therefore, the open boundary condition 
is a good approximation for calculating the transverse beam-beam force. 
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2.1 Green's Function 

Given a charge density Pc{x, y), which is normahzed to the total charge 

dxdypcix,y) = Ne, {2.1] 



where is the total number of particles, the electric potential y) satisfies the Poisson 
equation 

+ — U(x,i/) = -27rp,(x,y) (2.2) 



\dx'^ dy"^ 

with X and y being the transverse coordinates. The solution of the Poisson equation can be 
expressed as 

y) = j dx'dy'G{x -x', y- y')pc{,x\ y'), (2.3) 
where G is the Green's function which satisfies the equation 

^ +^ G{x -x\y- y') = -2n6{x - x')6{y - y'). (2.4) 



y dx"^ dy"^ 

In the case of open boundary condition, namely the boundary is far away so that its con- 
tribution to the potential can be ignored, one has the well-known explicit solution for the 
Green's function: 

G{x -x',y- yi) = In [(x - x'f + (y - y' f] . (2.5) 

This explicit solution can be used directly to compute the potential. The main problem of 
this approach is that it is slow to calculate the logarithm and the number of computations 
is proportional to the square of the number of macro particles Np. One can reduce Np by 
introducing a two-dimensional mesh to smooth out the charge distribution |Q . Or to further 
improve the computing speed, one can map the solution onto the space of spectrum by the 
Fast Fourier Transformation (FFT) and then calculate the potential |in . 



2.2 Reduce the Region of Mesh 

Another alternative approach is to solve the Poisson equation with a boundary condition 
0, because the region (20 pm x 450 pm for PEP-II) occupied by the beam is much smaller 
than the boundary defined by the beam pipe (2 cm radius) at the collision point. In order 
to achieve required resolution, a few mesh points per a of the beam are needed, otherwise 
the size of mesh is too large for numerical computation. 

However, it is unnecessary to cover the entire area with mesh inside the beam pipe since 
the area is mostly empty. We choose a smaller and finite area of the mesh, which is large 
enough to cover the whole beam, and by carefully selecting the potential on the boundary, 
we can obtain the accurate solution inside the boundary. 
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We denote by 0i the solution ( p.3| ) of the Poisson equation. Let 02 be the solution 
obtained by solving the Poisson equation in a two-dimensional area S with the potential 
prescribed on a closed one-dimensional L bounding the area S 



hix,y)= / dx'dy'G{x - x', y - y')pc{x\y'), 



(2.6) 



where (x, y) G L. By definition, we have (pi = 02 on the boundary L. Let U = 0i — 02 and 
use the first identity of Green's theorem |13[ in two dimensions 



/ 



/dU 
U-g^dl, 

L 



(2.7) 



where dl is a line element of L with a unit outward normal n. Since U = on L and V^?7 = 
inside L, we have 



y {yufdxdy = 0, 



(2.8) 



implying that U is a. constant inside L. We can set U = 0, which is consistent with the value 
on the boundary. Hence 0i = 02. The two solutions are identical. 



3 Field Solver 

We adopt the PIC technique to calculate the fields induced by the charge (and current) 
of the beams self-consistently. The charge distribution of a beam is represented by macro 
particles. These macro particles are treated as single electron or positron dynamically. In 
order to compute the field acting on the particles of the opposing beam, we first deposit 
their charges onto the gird points of a two-dimensional rectangular mesh. We denote by 
the horizontal distance between two adjacent grid points and by Hy the distance in vertical 
direction. 



3.1 Charge Assignment 



We choose the method of the triangular-shaped cloud |]15l as our scheme for the charge 
assignment onto the grid. On a two-dimensional grid, associated with each macro particle, 
nine nearest points are assigned with non- vanishing weights as illustrated in Fig. |l|. We use 
"0" to denote the first, "+" as the second, and "-" as the third nearest lines. 

The weights are quadratic polynomials of the fractional distance, = 6x/Hx, to the 
nearest line 



r 

4 ^ 

2(4^ 



(3.1) 
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Figure 1: Scheme of charge assignment. 



VI 2 

The coefficients are chosen such that the transition at the middle of the grid is continuous 
and smooth, and w'^ + + w~ = 1 which is required by the conservation of charge. In order 
to retain these properties, the weights of the two-dimensional grid are simply a product of 
two one-dimensional weights. For example, = w'^Wy or = w^w~. 

3.2 Poisson Solver 

It is crucial to solve the Poisson equation fast enough (within a second on a computer 
workstation) for the beam-beam simulation, because the radiation damping time is about 
5000 turns and several damping times are needed to reach an equilibrium distribution. For 
the reason of the computing speed, we follow Krishnagopal and choose the method of 
cyclic reduction and FFT WM. A five-point difference scheme is used to approximate the 



two-dimensional Laplacian operator 



-l,j + (Pi+l,j - ^(Pi,j , (Pi,j-1 + Vi,j + 1 - ^yij _ r,^^ /o o^ 

772 ^ 772 ~ -^T^PciJ, 1>5.ZJ 

X y 

where i and j are the horizontal and vertical indices that label the grid points on the mesh. 

Truncation errors are of the order of and Hy. It is worthwhile to mention that, if we 
use the same number of mesh points per a in both transverse directions in the case of beam 
aspect ratio 30:1, the truncation errors in the horizontal plane are dominant. To minimize 



the errors in our simulation, we select three times more mesh points per a in horizontal 
direction compared to the vertical one. 



3.3 Field 

^ — * 

The field E = — V0 is computed on the two dimensional grid, using a six-point difference 
scheme 

Exi,j = -^^^[{<Pi+i,o+i - <Pi-i,o+i) + - <i>i-i,j) + - ^i-ij-i)], (3.3) 

Eyi^j = - + - (t>i,3-l) + - (3.4) 

The field off the grid is computed with the same smoothing scheme used in the charge assign- 
ment to ensure the conservation of the momentum. The fields and Ey are interpolated 
between the grid points. They are calculated by using the weighted summation of the fields 
at the nine nearest points with exactly the same weights used as the charge is assigned. 

4 Track Particles 

The motion of a particle is described by its canonical coordinates 

^{x,P^,y,Py), (4.1) 
where Px and Py are particle momenta normalized by the design momentum po- 

4.1 One- Turn Map 

When synchrotron radiation is turned off, a matrix is used to describe the linear motion in 
the lattice 

Zn+l ^M-Zn, (4.2) 

where M is a 4 x 4 symplectic matrix which can be partitioned into blocks of 2 x 2 matrices 
when the linear coupling is ignored 

Here Mj., and My are 2x2 symplectic matrices. The matrix is expressed with the 
Cour ant- Snyder parameters P^, c^x, and at the collision point 

^ f cos(27rz/^) + aa.sin(27ri/^) sin(27rz/^) \ 
^ \ — 7a; sin(27rz/a;) cos(27ri/a;) — sin(27rz/^) i' 

where i/^ is the horizontal tune. A similar expression is applied in the vertical plane. 
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4.2 Damping and Synchrotron Radiation 



Following Hirata []I6|, we apply the radiation damping and quantum excitation in the nor- 
malized coordinates, since it is easily generalized to include the linear coupling. The motion 
of a particle in the normalized coordinate is described by a rotation matrix 

/ cos(27ri/^.) sin(27ri/^.) \ , . 

" 1^ - sin(27rz/,) cos(27ri/,) J ' ^ ' 

which is obtained by performing the similarity transformation 

= A-^ ■ M, ■ A^, (4.6) 



where 



A.= \ -J^ ^ \ .A-^=\ |. (4.7) 



When synchrotron radiation is switched on, we simply replace the rotation matrix R^ 
with following map in the normalized coordinates x and 



where rj^ and r]p^ are Gaussian random variables normalized to unity, is the damping 
time in unit of number of turns and ex is the equilibrium emittance. In the vertical plane, a 
similar map is applied. 

4.3 Beam-Beam Kick 

Assuming particles are ultra-relativistic and the collision is head-on, the kick on a particle 
by the opposing beam is given by the Lorenz force 



5P. = —etE., (4.9) 
SPy = -^Ey, (4.10) 

where and Ey are the horizontal and vertical components of the electric field evaluated 
at the position of the particle. They are computed with the Poisson solver as outlined in 
the previous section each time two slices of the beam pass each other. And the half of the 
transverse force is the magnetic force due the beam moving at the speed of light. The energy 
of the particle, Eq = cpQ, appearing in the denominator of the above expressions comes from 
the normalization of the canonical momenta P^ and Py and the use of the s-coordinate, 
s = ct, as the "time" variable. 



7 




Figure 2: The beam-beam kick by aflat Gaussian beam with aspect ratio 30:1 near X axis and 
Y axis. The dash- dotted curve is the case when (p = is assigned as the boundary condition. 
The long-dashed curve is the kick when inhomogeneous boundary condition is used. The 
short-dashed curve is the kick produced by the Erskine-Bassetti formula /fT/l/. 

A typical beam-beam kick experienced by a particle near the axis is shown in Fig. ^ 
with the PEP-II parameters, which are tabulated in the next section. As expected based on 
the derivation in section 2.2, the kick resulted from solving the Poisson equation with the 
inhomogeneous boundary condition agrees well with the analytic solution. In addition, the 
agreement demonstrates that the scheme of the charge deposition works well, the mesh is 
dense enough and the number of macro particles is large enough. 

The number of macro particles used to represent the distribution of the beam is 10240. 
The area of the mesh is 8(Ja;X24cry and there are 15 grid points per and 5 per ay. There 
are about 15 macro particles per cell within Sa of the beam. These parameters are chosen 
to minimize truncation errors and maximize resolution. The 256x256 mesh is also the 
maximum allowed by a computer workstation to complete a typical job within a reasonable 
time. 

The discrepancy between the solution with the homogeneous boundary condition, = 0, 
and the analytic one worsen as the beam aspect ratio becomes larger because the actual 
change of the potential on the horizontal boundaries becomes larger. 
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5 Simulation of PEP-II: Validation 



An object-oriented C++ class library has been written to simulate the beam-beam interac- 
tion using the method outlined in the previous sections. In the library, the beam and the 
Poisson solver are all independent objects that can be constructed by the user. For example, 
there is no limitation on how many beam objects are allowed in the simulation and the beams 
can have different parameters as an instance of the beam class. These features provide us 
with great flexibility to study various phenomena of the beam-beam interaction. 

We will carry out the simulation of beam-beam interaction with the current operating 
parameters of the PEP-II so that the results of the simulation can be compared with the 
known experimental observations. As a goal of this study, after a proper benchmarking of the 
code against the experiment, we shall be able to make predictions on parameter dependence 
and show how to improve the luminosity performance of the collider. 



5.1 PEP-II Operating Parameters 



Parameter 


Description 


LER(e+) 


HER( 


E (Gev) 


Beam energy 


3.1 


9.0 


P: (cm) 


Beta X at the IP 


50.0 


50.0 


P; (cm) 


Beta Y at the IP 


1.25 


1.25 


Tt (turn) 


Transverse damping time 


9740 


5014 


ex (nm-rad) 


Emittance X 


24.0 


48.0 


ey (nm-rad) 


Emittance Y 


1.50 


1.50 




X tune 


0.649 


0.569 


Uy 


Y tune 


0.564 


0.639 



Table 5.1: Parameters for the beam-beam simulation 



The parameters used in the simulation are tabulated in Tab. |5.1| . The vertical /?* is 



lowered to 1.25cm [0 from the design value 1.5cm |TT[. The horizontal emittance 24nm-rad 
in the Low Energy Ring (LER) is half of the design value 48nm-rad because the wiggler is 
turned off to increase the luminosity. The damping time, 9740 turns, in the LER is a factor 
of two larger than the one in the High Energy Ring (HER) because of the change of the 
wigglers made during the construction of the machine. The degradation of luminosity from 
the increase of the damping time was found then to be about 10% based on the beam-beam 
simulation. The tunes are split and are determined experimentally to optimize the peak 
luminosity. 



5.2 Procedure of simulation 

The distribution of the beam is represented as a collection of macro particles that are dy- 
namically tracked. The procedure to obtain equilibrium distributions of the two colliding 
beams is as follows 
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• initialize the four-dimensional Gaussian distribution according to the parameters of the 
lattice at the colhsion point and the emittance of the beam. Distributions of two beams are 
independent and different. 

• iterate a loop with three damping times 

• propagate each beam through corresponding lattice using one-turn map with syn- 
chrotron radiation. 

• cast the particle distributions onto the grid as the charge distribution with weighting 
and smoothing. 

• solve for the potential on the grid with the Poisson solver. 

• compute the field on the grid. 

• calculate the beam-beam kick to the particles of the other beam with the field at the 
position of the particles. The field off the grid is interpolated with the same weighting and 
smoothing used in the charge deposition. 

• save data such as beam size, beam centroid and luminosity. 

• end of the loop. 

• save the final distributions. 

We vary the beam intensity with a fixed beam current ratio: /_|.:/_ = 2:1, which is close 
to the ratio for the PEP-11 operation. At each beam current, we compute the equilibrium 
distributions. 

5.3 Beam-Beam Limit 

Given equilibrium distributions that are close enough to the Gaussian, we can introduce the 
beam-beam parameters 



= (5.1) 



where rg is the classical electron radius, 7 is the energy of the beam in unit of the rest 
energy, and is total number of the charge in the bunch. Here the superscript "-f" denotes 
quantities corresponding to the positron and "-" quantities corresponding to the electron. 

The results of the simulation are shown in Fig. |^. The beam-beam tune shifts for the 
electron beam are low because of the large beam-beam blowup of the positron beam. At this 
operating point, the positron is the weaker beam. When /_|_ = 1200mA and /_ = 600mA, 
which is the near the maximum allowed currents when the beams are in collision, the positron 
beam sizes are = 260/im and = 7fim. 

5.4 Luminosity 

Given the two beam distributions, and p~ , the luminosity can be written as 
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Figure 3: The beam-beam tune shifts as a function of beam currents. Number of bunches, 
Ub = 554, is used for the total beam currents. The revolution frequency /o = 136.312 kH^ 



L = nbfoN^N J J p^{x,y)p {x,y)dxdy, (5.2) 

— oo — oo 

where Ub is the number of the colliding bunches, /o is the revolution frequency, and A^"*", 
are the number of charges in each position and electron bunch, respectively. Since the 
distribution p is normalized to unity 

J dxdyp{x,y) = 1 (5.3) 

and proportional to the charge density pc, we evaluate the overlapping integral by a summa- 
tion over p^p~ on the mesh. If we assume the distributions are Gaussian, the overlapping 
integral can be carried out 
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where S^; = y(T+ + cr~ and Tiy = Ja+ + • Two methods agree within a few percents. 
The mesh method gives a higher luminosity than the Gaussian one. We always use the mesh 
method, since it can be applied to broad classes of distribution. 
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Figure 4: Luminosity as a function of the beam current. The labels are the number of the 
colliding bunches. 

Figure ^ shows the luminosity of the beams with 415 colliding bunches, which are spaced 
with every 8 RF buckets and 10% of the gap. The luminosity is beam-beam limited. It also 
shows that the optimum number of bunches is between 544 and 665 and the luminosity 
is about 2.3xl0^^cm~^s~^ given J+ = 1200mA. These results quantitatively agree with 
the experimental observations in the routine operation of the PEP-II. For example, the 
peak luminosity of the PEP-II is 1.95x lO^^cm'^s"^ with /+ = 1170mA, /_ = 700mA, 
and Uf, = 665 during the period of June, 2000. The fact that the luminosity value in the 
simulation is higher than the observation could be explained by the hour-glass effect which 
is ignored in the simulation. 



12 



For a fixed number of bunclies, say 554, tlie simulation sfiows a maximum luminosity, 
which is also seen daily in the control room of the PEP-II. From the simulation, we see that 
the reason for the peaked luminosity is the rapid growth of cfy once the peak current is 
passed. 

In addition, the simulation predicts that we can reach the design luminosity 3x 10'^'^cm~^s~^ 
by running 829 bunches at the beam current of /+ = 1600mA and /_ = 800mA. This predic- 
tion has not been realized yet at this time. Currently, the total positron current is probably 



limited below 1200mA by the electron-cloud instability |19|. Once this limitation is removed 



we expect to reach the design luminosity with 829 bunches. 

There is no particle loss outside the area (8cr2^x24(Ty) covered by the mesh in the first 15 
data points. Beyond the 15th points, particle loss is almost about 1%. 



5.5 Damping Time 

Historically, the damping time is typically not considered to be an important parameter for 
the beam-beam effects. So we make an attempt to reduce the damping time artificially for 
the LER to speed up the computation. The result is shown in Fig. ^ 

The only difference of the parameters used in two simulations is the damping time in the 
LER, which is indicated as the labels in the figure. Indeed, at the low current, the difference 
of the luminosity is rather small, which is consistent with the simulation performed when the 
change of the wiggler was made. But the difference grows larger, as the current increases. 
At the peak luminosity for the PEP-II operation, = 1200mA, the difference is about 40%, 
which is significant. 

This result shows for the first time that the damping time is a rather important parameter 
for the computation of the peak luminosity at high beam currents. Secondly, it points a way 
to improve the peak luminosity of the PEP-II without the increase of the beam currents, 
namely to install another wiggler in the LER to reduce the damping time to the original 
design value. 



6 Discussion 

We have developed a hybrid method of solving the potential with an open boundary by 
using Green's function to fix the potential on a finite boundary and then to solve the Poisson 
equation for the potential inside the boundary. The method is applied to the simulation of 
strong-strong interaction of beam-beam effects in PEP-II. The preliminary results of this 
simulation show a very good quantitative agreement with the experimental observations. 
Given the simplicity of the two-dimensional model used, the achievement is surprising and 
remarkable. We have demonstrated that the present code has a highly reliable predictive 
capability of realistic beam-beam interaction. To further benchmark the code, we need to 
extend the simulation to include the finite length of the bunch and compare the simulation 
results directly to the controlled experiments. 

This method is quite general. It can be applied to the problem of space charge in three 
dimensions. It can also be used in the beam-beam interaction of a linear collider. Finally, 
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Figure 5: Luminosity affected by the damping time with 554 bunches. 



it can be applied to any boundary condition to reduce the region of the mesh if Green's 
function is known. 
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